*!version 1.4.1  01/07/24

mata
mata drop comb_n_local()
void comb_n_local() {
	//# setup 
		// setup timer
		timer_clear()
		timer_on(1)	
		
		// setup parties, party-ids, and portfolios
		mat_portfolios_parties = (st_data(1, "party_portfolios*")\st_data(1, "party_id*")\(1..cols(st_data(1, "party_portfolios*"))))
		vec_select = colmissing(mat_portfolios_parties):==0
		mat_portfolios_parties = select(mat_portfolios_parties, vec_select)
		sca_parties_num = cols(mat_portfolios_parties)
		
		// setup information on coalition parties
		vec_party_coalition = select((st_data(1, "party_coalition*")),vec_select)
		
		// setup information on mayor party
		vec_party_mayor = select((st_data(1, "party_mayor*")),vec_select)
		
		// setup seat-share
		vec_seat_share = (st_data(1, "party_seat_share*"))
		vec_seat_share = editmissing(select(vec_seat_share, vec_select),0)
		vec_seat_share = vec_seat_share:/sum(vec_seat_share)
		
		// setup utilities based on acutal actual allocation
		vec_utilities = select((st_data(1, "party_portfolio_share*")),vec_select)			
		
		// setup party-specific saliencies and their sum
		mat_party_saliencies = (st_data(1, "portfolio_salience*_party*"))
		mat_party_saliencies = rowshape(mat_party_saliencies,cols((st_data(1, "portfolio_salience*_party1"))))[.,1..sca_parties_num]'
		vec_select = colmissing(mat_party_saliencies):==0
		mat_party_saliencies = select(mat_party_saliencies, vec_select)
		mat_party_saliencies_total = rowsum(mat_party_saliencies)
		mat_party_saliencies_total[selectindex(mat_party_saliencies_total:==0)] = J(length(selectindex(mat_party_saliencies_total:==0)), 1, 1234)
		mat_party_saliencies = mat_party_saliencies:/mat_party_saliencies_total
		sca_portfolios_num = cols(mat_party_saliencies)
		
		// setup actual allocation
		vec_allocation_tmp = select(st_data(1, "portfolio_party_id*"), vec_select)	
		vec_allocation = (vec_allocation_tmp:==mat_portfolios_parties[2,1]):*1
		for (p=2; p<=cols(mat_portfolios_parties); p++) {
			vec_allocation = vec_allocation :+ (vec_allocation_tmp:==mat_portfolios_parties[2,p]):*p
		}
		
		// setup ID
		sca_id = (st_data(1, "id"))
		
		// setup filenames
		sca_filename = strofreal(sca_parties_num)+"parties_"+strofreal(sca_portfolios_num)+"portfolios"
		sca_filename
		
		// sca_portfolios_num = cols(vec_allocation)
		mat_portfolios_parties
		mat_party_saliencies
	//# setup end

		
	//# generate all potential allocations
		// generate all compositions, i.e. n as the sum of a sequence of (strictly) positive integers. 
		mat_compositions = mat_compositions_select = (mm_compositions(sca_portfolios_num,sca_parties_num))'
		mat_partitions = mm_partitions(sca_portfolios_num,sca_parties_num)
		sca_compositions = rows(mat_compositions)
		vec_out = rowsum(mat_compositions:==mat_portfolios_parties[1,]):==sca_parties_num
		// -> mat_compositions contains all possible compositions
		// draw random subsample if all possible compositions are too many
		if (sca_compositions >= 30000) {
			sca_subset = 2000/rows(mat_compositions)
			
			vec_actual_alloc = select(mat_compositions, vec_out)
			rseed(1234)
			vec_select = runiform(rows(mat_compositions), 1) :< sca_subset
			mat_subset_alloc = select(mat_compositions, vec_select)

			mat_compositions = mat_compositions_select = vec_actual_alloc \ mat_subset_alloc
			sca_compositions = rows(mat_compositions)
		}
	//# generate all potential allocations end
	
	//# estimate office payoffs for each potential allocation
		mat_office_payoffs = mat_compositions:*(1/sca_portfolios_num)
		// estimate proportionality for each potential allocation 
		vec_prop = 1:-(rowsum(abs(mat_office_payoffs:-vec_seat_share)))
		// estimate equitability for each potential allocation 
		vec_equi = 1:-(rowsum(abs(mat_office_payoffs:-J(1,sca_parties_num,(1/sca_parties_num)))))		
		// indicator variable for actual outcome
		vec_out = rowsum(mat_compositions:==mat_portfolios_parties[1,]):==sca_parties_num
		// indicator variable for potential allocation only coalition parties 
		vec_coal = rowsum((mat_compositions:>0):==vec_party_coalition):==sca_parties_num
		// indicator variable for allocation including mayor party
		vec_mayor = select(mat_compositions:>0, vec_party_mayor)
		if (rowsum(vec_party_mayor) == 0) vec_mayor = J(sca_compositions,1,0)
		
		mat_res = J(sca_compositions,1,sca_id),vec_out,vec_prop,vec_equi,vec_coal,vec_mayor
	//# estimate office payoffs for each potential allocation end
	
	//# estimate policy payoffs for each potential allocation
		if (sca_compositions < 30000) {
		p_vec_compositions = J(sca_compositions,1,NULL)
		sca_stored = 0
		// -> p_vec_compositions collects results for all compositions
		
		// check if local file exists
		if (fileexists(sca_filename)) {
			fh_get = fopen(sca_filename, "r")
			p_vec_compositions = fgetmatrix(fh_get) 
			fclose(fh_get)
			sca_stored = 1
		}				
		
		// loop through compositions	
		printf("%7.0f total compositions \n", sca_compositions) 
		if (sca_stored==1) printf("using stored matrices \n")
		if (sca_stored==0) printf("generating new matrices \n")
		
		for (composition=1; composition<=rows(mat_compositions); ++composition) {		
			
			printf("composition: %7.0f \n", composition) 
						
			mat_portfolios_parties = select(mat_compositions[composition,],mat_compositions[composition,]:>0)
			vec_parties = mm_which(mat_compositions[composition,]:>0)
			sca_parties_num = cols(mat_portfolios_parties)
		
			sca_comb = 1
			vec_comb_num = J(1,cols(mat_portfolios_parties)-1,.)
			vec_pos_max = J(1,cols(mat_portfolios_parties)-1,0)
			if (sca_parties_num == 2) {
				sca_comb = comb(rowsum(mat_portfolios_parties[1,1..cols(mat_portfolios_parties)]),mat_portfolios_parties[1,1])
				vec_comb_num = sca_comb
			}
			else {
				for (p=1; p<=sca_parties_num-1; p++) {
					vec_comb_num[.,p] = comb(rowsum(mat_portfolios_parties[1,p..cols(mat_portfolios_parties)]),mat_portfolios_parties[1,p])
					sca_comb = sca_comb * vec_comb_num[.,p]
					vec_pos_max[.,p..cols(vec_pos_max)] = vec_pos_max[.,p..cols(vec_pos_max)]:+J(1,(cols(vec_pos_max)-p+1),mat_portfolios_parties[1,p])
				}
			}
			
			if (sca_comb < 1000000) {
			
			mat_res_sal = J(sca_comb,sca_portfolios_num,0)
			mat_alloc = J(sca_comb,sca_portfolios_num,0)
			
			if (sca_parties_num > 2) {
				vec_mdelta = J(1,cols(mat_portfolios_parties)-2,.)
				vec_mdelta[1,1] = sca_comb/vec_comb_num[1,1]
				for (p=2; p<=sca_parties_num-2; p++) {
					vec_mdelta[1,p] = vec_mdelta[1,p-1]/vec_comb_num[1,p]
				}
			}
			
			if (sca_stored==0) {
				// populate results matrix and loop over parties: generate n over k combinations
				for (p=1; p<=sca_parties_num-1; p++) {
					
					printf("party: %7.0f \n", p) 
					
					k = mat_portfolios_parties[1,p]
					n = rowsum(mat_portfolios_parties[1,p..sca_parties_num]) 
					sca_comb_num = comb(n,k)
					
					for (i=1; i<=sca_comb_num; i++) {
						X = i - 1
						nfound = 0
						nmax = n - 1
						vec_comb = J(1,k,.)
						while(nfound < k) {
							kact = k - nfound
							nact = nmax
							while (comb(nact,kact) != . & comb(nact,kact) > X) nact = nact - 1
							nfound = nfound + 1
							vec_comb[1,nfound] = nact
							X = X - comb(nact,kact)
							nmax = nact - 1
						}
						vec_comb = vec_comb:+1
						if (p==1) {
							sca_mmin = (i-1)*(sca_comb/vec_comb_num[.,p])+1
							sca_mmax = i*(sca_comb/vec_comb_num[.,p])
							mat_res_sal[sca_mmin..sca_mmax,vec_comb] = J((sca_comb/vec_comb_num[.,p]),mat_portfolios_parties[1,p],vec_parties[p])
						}
						if (sca_parties_num > 2) {
							if (p>1 & p<(sca_parties_num-1)) {
								sca_mrange = vec_comb_num[1,p+1]
								for(m=p+2;m<=cols(vec_comb_num);m++) {
									sca_mrange = sca_mrange*vec_comb_num[1,m]
								}
								sca_mdelta = vec_mdelta[1,p-1]
								sca_mmin = ((i-1)*sca_mrange)+1
								sca_mmax = i*sca_mrange
								vec_mrange = (range(0, sca_comb-sca_mdelta, sca_mdelta)):+1
								vec_mrange = J(rows(vec_mrange),1,(sca_mmin,sca_mmax)):+vec_mrange:-1
								for (j=1;j<=rows(vec_mrange);j++) {
									mat_alloc[vec_mrange[j,1]..vec_mrange[j,2],(vec_pos_max[.,p-1])+1..vec_pos_max[.,p]] = J((vec_mrange[j,2]-vec_mrange[j,1]+1),1,vec_comb)
									mat_avail = mm_which(mat_res_sal[vec_mrange[j,1],.]:==0)
									mat_res_sal[vec_mrange[j,1]..vec_mrange[j,2],mat_avail[.,vec_comb]] = J((vec_mrange[j,2]-vec_mrange[j,1]+1),mat_portfolios_parties[1,p],vec_parties[p])
								}	
							}
							if (p==(sca_parties_num-1)) {
								sca_mrange = range(i,sca_comb-vec_mdelta[1,p-1]+i,vec_mdelta[1,p-1])
								for (jj=1;jj<=rows(sca_mrange);jj++) {
									mat_avail = mm_which(mat_res_sal[sca_mrange[jj,1],.]:==0)
									mat_res_sal[sca_mrange[jj,.],mat_avail[.,vec_comb]] = J(1,mat_portfolios_parties[1,p],vec_parties[p])
								}
							}
						}
					}				
				}
				
				mat_res_sal = mat_res_sal:+((mat_res_sal:==0):*vec_parties[sca_parties_num])
				p_vec_compositions[composition] = &(mat_res_sal:+((mat_res_sal:==0):*sca_parties_num))				
			} 
			
			else if (sca_stored==1) {
				mat_res_sal = *(p_vec_compositions)[composition]
			}
			
			mat_port_party_sum = J(rows(mat_res_sal), rows(mat_party_saliencies), 0)
			for(party=1; party<=cols(vec_parties); party++) {
					mat_port_party_sum[.,vec_parties[party]] = rowsum((mat_res_sal:==vec_parties[party]) :* J(rows(mat_res_sal),1,(mat_party_saliencies[vec_parties[party],.])))
			}
						
			if(composition == 1) {
				mat_res_all = J(rows(mat_res_sal),1,mat_res[composition,]),(1:-rowsum(abs(mat_port_party_sum:-vec_seat_share))),(rowsum(mat_res_sal:==vec_allocation):==cols(mat_res_sal))
			}
			else {
				mat_res_all_add = J(rows(mat_res_sal),1,mat_res[composition,]),(1:-rowsum(abs(mat_port_party_sum:-vec_seat_share))),(rowsum(mat_res_sal:==vec_allocation):==cols(mat_res_sal))
				mat_res_all = mat_res_all \ mat_res_all_add
			}
			
			}
			
		}
	//# estimate policy payoffs for each potential allocation end
	
		stata("clear")
		st_addobs(rows(mat_res_all))
		idx = st_addvar("double", ("id","outcome_quant","proportionality","equitability","coalition","mayor_party","qual_prop","outcome_qual"))
		st_store(.,1..cols(mat_res_all), mat_res_all) 
		}
		else {
			stata("clear")
			st_addobs(1)
			idx = st_addvar("double", ("id","outcome_quant","proportionality","equitability","coalition","mayor_party","qual_prop","outcome_qual"))
			st_store(.,1..8, (sca_id,J(1,7,0))) 		
		}
		// stop timer		
		timer_off(1)
		timer_value(1)[1,1]
		
}
end	
